function [CE,P_star,P_star_0] = Get_CE_TDF_optim_fun(P_choices,MF_fee_rate,u,T,g_I,R_S_GHQ,weight,r,A_grid,I_0)

%% Set model specifications & parameters
P_num = length(P_choices);


%% State space
A_num = length(A_grid); A_grid_max = max(A_grid);   % grid of potential QRP account values A_t.


%% Matrizes to store results from dynamic optimization
EU_star = zeros(A_num,T);           % to store expected utility of being in current state (and time)
P_star = zeros(A_num,T-1);          % to store optimal equity exposure rate


%% Terminal utility
EU_star(:,T) = u(max(A_grid,0.01));


%% Recursion
for t = (T-1):-1:1
    EU_tplus1 = EU_star(:,t+1);     % value function at end of period.
    I_t = I_0 * (1+g_I)^t;  % PH's time-t contribution to QRP (no contribution at retirement time T).
    for a = 1:A_num
        A_t = A_grid(a);

        % Analyze all potential equity exposure rates:
        EU_P = u(0.001) * ones(P_num,1);
        for p = 1:P_num
            P = P_choices(p);
            R_A = max( P * R_S_GHQ + (1-P) * (exp(r)-1) , -1 );             % GHQ nodes for return credited to QRP account.
            A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-MF_fee_rate) , A_grid_max );         % GHQ nodes for QRP account value at year-end.
            U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);              % GHQ nodes for utility at year-end.
            EU_P(p) = weight' * U_tplus1;                               % expected utility of account values at end of period.
        end

        % find value of 'P' that maximizes expected utility in current state (& time):
        EU_max_P = max(EU_P);
        pp = find(EU_P == EU_max_P);
        P_star(a,t) = P_choices(pp(1));
        EU_star(a,t) = EU_max_P;
    end
end


%% Time 0

EU_tplus1 = EU_star(:,1);     % value function at end of period.
I_t = I_0;      % PH's time-t contribution to QRP (no contribution at retirement time T).
A_t = 0;

% Analyze all potential equity exposure rates:
EU_P = u(0.001) * ones(P_num,1);
for p = 1:P_num
    P = P_choices(p);
    R_A = max( P * R_S_GHQ + (1-P) * (exp(r)-1) , -1 );             % GHQ nodes for return credited to QRP account.
    A_tplus1 = min( (A_t + I_t)*(1+R_A)*exp(-MF_fee_rate) , A_grid_max );         % GHQ nodes for QRP account value at year-end.
    U_tplus1 = interp1(A_grid,EU_tplus1,A_tplus1);              % GHQ nodes for utility at year-end.
    EU_P(p) = weight' * U_tplus1;                               % expected utility of account values at end of period.
end

% find value of 'P' that maximizes expected utility in current state (& time):
EU_max_P = max(EU_P);
pp = find(EU_P == EU_max_P);
P_star_0 = P_choices(pp(1));
EU_star_0 = EU_max_P;


%% Convert utility to Certainty Equivalent (CE)
options = optimset('Display','none');
CE = fsolve(@(ce)  1e+3 * (u(ce)/EU_star_0-1), 200,options);    % certainty equivalent to potential values of A_T.